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Abstract 

We present a simple parameterization of the running coupling constant av{q), 
defined via the static potential, that interpolates between 2-loop QCD in the UV 
and the string prediction V{r) = o'r — ^ in the IR. Besides the usual A- 
parameter and the string tension a, av{q) depends on one dimensionless parameter, 
determining how fast the crossover from UV to IR behavior occurs (in principle we 
know how to take into account any number of loops by adding more parameters). 
Using a new ansatz for the lattice potential in terms of the continuum av{q), we 
can fit quenched and unquenched Monte Carlo results for the potential down to 
one lattice spacing, and at the same time extract av{q) to high precision. We 
compare our ansatz with 1-loop results for the lattice potential, and use ay(g) 
from our fits to quantitatively check the accuracy of 2-loop evolution, compare 
with the Lepage-Mackenzie estimate of the coupling extracted from the plaquette, 
and determine Sommer's scale tq much more accurately than previously possible. 
For pure SU(3) we find that ay(g) scales on the percent level for /3 > 6. 



PACS numbers: 12.38.Aw, ll.15.Ha, 12.38.Gc 



1 Introduction 



In a situation where quarks are non-relativistic their dynamics is much slower 
than that of the massless gluons. In the limit of infinitely heavy quarks, or, equiv- 
alently, static sources, their interaction can be described by a local, instantaneous 
potential V{r). This static potential will for the purposes of this paper always 
refer to the potential between a quark and an anti-quark in the fundamental rep- 
resentation of SU(A^). 

For very short distances, V{r) is Coulomb like, with a running charge that is 
known to two loops in terms of the (unknown) A-parameter. The strong coupling 
expansion of lattice gauge theory suggests a string picture, where a flux tube 
connects two widely separated quarks, in accord with the idea of confinement [2] 
(we here ignore the string breaking in the presence of dynamical fermions). In 
addition to a linear term involving the string tension a, the potential has a universal 
large r correction [3] due to the zero point energy of the transverse fiuctuations of 
the flux tube, so that in the IR 



up to a non-universal constant. The "IR charge" Bir = 7i{D — 2)/24 in D space- 
time dimensions, assuming the transverse fluctuations are described by the simplest 
bosonic string theory.^ 

At present the only way to obtain quantitative information on V{r) is through 
the Monte Carlo simulation of Wilson loops. These are nowadays very precise; 
the lattice potential can be measured with a relative accuracy approaching 10~^ at 
short distances and 10^^ at long distances. The lattices used are typically of size 
32^ for pure gauge theory and 16^ in the presence of dynamical quarks. The spacing 
of the flnest lattices available for SU(3) is about 0.04 fm, corresponding to 5 GeV, 
just where 2-loop perturbation theory is expected to become good. However, at 
these distances the potential suffers from large lattice artifacts, which obscure the 
view of perturbative QCD. These lattice artifacts are on the order of 10% at one 
lattice spacing if one uses the standard (unimproved) Wilson action for the gluons, 
and fall roughly like at larger distances. One therefore has to go to many 
lattice spacings before rotational invariance is restored within errors. 

In early studies, the lattice potential, at least restricted to on-axis points, was 
fitted to an ansatz of the form —e/r + ar + Vo, i.e. it was assumed to behave as in 

^The situation is more subtle in the presence of a non-perturbative gluon condensate. But 
it is now clear that the presence of a gluon condensate is not necessarily incompatible with the 
existence of a local potential. See [1] for a review and references. 

^Nowadays [4] this universal correction is easily recognized as the Casimir energy of the 
(bosonic) conformal field theory in a finite geometry with free boundary conditions that describes 
the transverse fluctuations of the flux tube; D — 2 being the central charge of this theory. 



V{r) = ar — — + . . . , 
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the IR, eq. (1). Not just does this ansatz completely ignore the lattice artifacts, 
it strictly speaking contradicts QCD at short distances. The number e has to be 
interpreted as some kind of effective average coupling. It is not even clear that e is 
an effective short-distance quantity, since it usually comes out rather close to the 
IR charge e^^ = 7r/f2 in eq. (f). 

To partially accommodate lattice effects, Michael [5] has suggested to add a 
term to this ansatz, proportional to the difference between the lattice and the 
continuum Coulomb potential. This takes into account at least the tree-level lattice 
artifacts. Also, since recent Monte Carlo data [5, 6, 7] on fine lattices arc accurate 
enough to see that e really is a running charge, Michael [5] proposed to take this 
into account in a phenomenological way by replacing e with e — //r, with some 
new parameter /. If necessary, we will distinguish these kind of ansatze from the 
naive Coulomb -|- linear one, by calling them "modified Coulomb -|- linear ansatze" 
(or fits). 

With this setup one can include smaller r- values in the fit, but the first few 
still have to be left out — and more and more as the data get more precise and the 
lattice spacings smaller. Furthermore, the short-distance parameters are not very 
stable with respect to varying the fit range, in particular on fine lattices. This 
is not surprising, since asymptotic freedom is not properly incorporated in this 
ansatz; the short-distance parameters are all effective parameters. Therefore, the 
only physical result one obtains from these fits is, at best, an (improved) estimate 
of the string tension a; no information about short-distance QCD, e.g. the A- 
parameter, is gained. 

To learn about A another idea of Michael [5] is now widely used. Namely, one 
defines a running coupling apir) in the "force scheme" by 

rV(r) = C^a^(r) . (2) 

Here Cp = (iV^ — 1)/2N is the Casimir invariant of the fundamental representation 
of SU(A^). One then fits numerical derivatives of the potential data to a 2-loop 
formula for apir). Lattice artifacts, errors from taking numerical derivatives, and 
the fact that 2-loop evolution is only on the verge of becoming reliable, lead to 
quite large errors in A. Furthermore, it is hard to estimate the systematic errors 
of this A-determination; one might simply be determining an effective A, that is 
still quite far from the true A. This could be true even if some consistency checks 
involving different lattice spacings, e.g. scaling of quantities involving A, work quite 
well. In fact, our results indicate that this is exactly what happens. 

We here describe a new method that allows for a unified fit of potential data 
at all distances by incorporating lattice artifacts and the running of the coupling 
in a more fundamental way. The general idea is to express to the lattice potential 
in terms of a running coupling. More precisely there are three ingredients: 
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1. As advocated in refs. [8, 9], for example, use the continuum static potential 
to define a physical running coupling CKv{q) via 

^ . (3) 

This scheme for the coupling will be referred to as the V scheme. With this 
definition the effect of multi-gluon exchange is absorbed into av{q)- Fourier 
transforming eq. (3), the continuum potential is by definition given by 

p d^g _iqr 47ray(g) 2Cf , singr 

(4) 

2. Globally parameterize av{q) to take into account 2-loop QCD in the UV in 
terms of A, the string prediction of eq. (1) in the IR in terms of a, and to 
have another parameter determining how fast the crossover from the UV to 
the IR behavior occurs. [Actually we will see that by introducing further 
parameters we can take into account any number of loops in the UV.] 

3. Finally, make the ansatz that in terms of the continuum av{q) the potential 
on an infinite lattice of spacing a can be accurately represented as, 

Va^{r) = Vo - Cf -— ^ e '^i'^ , g^ = - sm — , 5 

J-TT/a {Zny q^ a 2 

except presumably for r = 0, where the whole underlying picture of gluon 
exchange breaks down (which is why we have to add the constant Vq after 
making the integral convergent by subtracting its value at r = 0).^ 

The three parameters in av{q) and the constant Vq are then fitted by matching 
eq. (5) to the Monte Carlo data for the potential. 

To be sure, the above equation for the lattice potential in terms of the contin- 
uum coupling must be regarded as an ansatz. The motivation for this ansatz is 
that it appears to be the most natural expression that (a) gives the exact tree level 
lattice potential, that is, the lattice Coulomb potential (corresponding to the case 
av{q) = const), and (b) has the right continuum limit. The fact that eq. (5) will 
allow us to obtain excellent fits of the lattice potential of a wide variety of theories 
down to one lattice spacing, provides indirect empirical evidence for its accuracy. 
Additional cnideiiee will he discussed later. 



^The "lattice momentum" q in eq. (5) is appropriate for the usual Wilson action. If an 
improved action is used, q should be replaced by the corresponding "improved momentum", 
yielding an improved potential Va* (r) , etc. 
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The outline of this paper is as follows. In sect. 2 we provide the details of 
the parameterization of av{q) satisfying the required UV and IR behavior. In 
particular, we prove an elementary theorem showing how the running coupling 
can be written in terms of "iterated logarithms" to all orders of the loop expansion 
of the /9-function. This theorem then suggests a natural way of avoiding the 
perturbative Landau pole, motivating our ansatz for the running coupling. 

In sect. 3 we discuss the efficient evaluation of the three-dimensional integral in 
eq. (5), and details of our fitting and error analysis procedure. In sect. 4 we apply 
our results to Monte Carlo data for the lattice potential of various gauge theories; 
with gauge group SU(3) or SU(2), with or without dynamical fermions (Wilson 
and staggered). Using eq. (5) leads to excellent fits of the data down to one lattice 
spacing. To check for systematic errors in the running coupling av{q) extracted 
from these fits, we perform fits with an ay(g) that partially incorporates 3-loop 
effects. We compare our results with the Lepage-Mackenzie estimate [10] of av{q) 
in the UV, in particular also on very fine lattices. We then compare the 1-loop 
expansion of eq. (5) with the 1-loop results of [11] for the on-axis lattice potential. 
Curiously enough, our 1-loop expansion agrees much better with the Monte Carlo 
data than the results of [11]. These comparisons provide good evidence for the 
accuracy of representing the lattice potential as in eq. (5). Finally, we check how 
accurate 2- loop evolution is compared to our non-perturbative av{q)- 

Sect. 5 is devoted to scahng. We use the intermediate scale Tq introduced by 
Sommer [12] to determine the ratios of lattice spacings of various theories much 
more accurately than previously possible. For pure SU(3) we find that ay (g) scales 
on the percent level for /? > 6.0. In sect. 6 we compare our approach with other 
schemes to define and extract a running coupling from Monte Carlo simulations. 
In particular, we provide some details about the force scheme and its close relation 
to the V scheme. We conclude in sect. 7 by summarizing our findings and outlining 
various applications and extensions. 

In [13] we presented a brief summary of our results. Due to the fact that our 
error analysis is somewhat involved, some of the results given in [13] were slightly 
preliminary (also, in one case we only had preliminary potential data). 

2 Parameterizing the Running Coupling 

2.1 The Ultraviolet Regime 

We will first describe some general results concerning the integration of the 
function 

l3{a) = q^dq2a{q) = —(Soa^ — (3ia^ — (320.^ — ... , (6) 
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of the coupling a. — g'^/47r in a generic scheme. Recall that the 1- and 2-loop 
coefficients are scheme independent and equal to [14] 

4vr/?o = yA^-^n,, (4;rr = | AT^ - (^^TV + ^) (7) 

in QCD with gauge group SU(A^) and flavors of masslcss fcrmions. The higher 
coefficients depend on the scheme (none are presently known in the V scheme). In 
terms of Q;(g) = /3oQ;(g) the integration of the /3-function immediately leads to 



- = i + 61n — — -/ dxU— + ^-- — =i + 61n- 



where /9(q;) = -q;^(1 + {5ia + (320l^ + ■■■)> A = Pi/Po , and 

t = ln(|^) , 6 = A , h = P!-P2 , h = PA-hP!-A ■ (9) 

We have followed the usual convention of defining the integration constant, the 
A-parameter, so that there is no constant term on the rhs of (8). Note that as long 
as the /3-function has a power series expansion in a, there is only one logarithmic 
term on the rhs of (8), coming from the first two loops. Solving (8) iteratively 
gives 



1 , n 1 

- = t + b\nt + - 

a t 



1 



1 



bnnt + b2 +^ --bnnU + {b-^-bb2)lnt + bb2 + b3 +.... (10) 
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As an aside we would like to point out that the full 2-loop expression for 1/a 
involves not just the leading terms t+blnt usually quoted, but also the subleading 
2-loop term b^lnt/t. 

Recall that the above expansion has a Landau pole, i.e. the perturbative a{q) 
diverges for some q > A as q approaches A from above. One of the goals of the 
ansatz we will later use for av{q) is of course to avoid this unphysical Landau pole. 
To motivate our ansatz we now describe a different way of solving (8) for a{q). 
Namely, consider the following recursively defined function 

R(t) = R^''\t) , R^''\t) = b\n[ak{t + R^''+^\t))] , k^0,l,'2,.... (11) 

We now prove that by suitably choosing the coefficients > the running coupling 
can be written as 

- = t + R{t) . (12) 

a 

This is equivalent to showing that R{t) = t — ^ can be rewritten in terms of a 
in the form R — bla^ + b2a + 630;^ -|- The proof is simple: By substituting 
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i = 4 - in the definition of R^^'^ we have 

a 

= 6 In (^) + 6 In [l + - R^'^^)] 



(13) 



All we have to do is to iterate these equations, starting with R^^'^ = 61n(^). By 
induction one therefore proves 



n 



(14) 



n=l 



where R!^^ — lnajfc+„ + terms involving only Inoj with j < n + k. This completes 
the proof. 

Specifically, = ln(ajk+i/ao), -^2''^ = ln(afe+2/ai) - | ln^(ajt+i/ao), etc. So 
the first few are determined by 

Oo = 1 , ft^lnoi = 62, 6^(lna2 — Inoi — ^ In^ ai) — 63. (15) 

Note that to obtain the 1/d expansion to n loops one can truncate the iterative 
definition of R{i) by setting i?^"""^) = feint. For example, to get two loops we can 
simply set R{t) = 61n(i + feint), cf. (10). 

Our iterated form of R{i) still has a Landau pole, but now a simple way to 
avoid it suggests itself by the replacements: 



1 



Q;(g) 



In 



1 + e* e^'°^W 



and 



In" 



(16) 



(17) 



with some constants > 1. [The constant 1 in eq. (16) is chosen to get the right 
IR behavior for ay{q)^ see below. In a scheme where the couphng "freezes out" in 
the IR some constant > 1 should be chosen.] Note that the Cjt are non-perturbative 
parameters. Obviously the new «(g) has no Landau pole. 
To two loops we can now write as ansatz for Q;y(g) 



/?oa;y(g) 



In 



1 + e* ln'[co + e*A(t)] , h = (3J(5l , t = ln(gVA 



where we know that InA(i) = R^^\t) — bint in the UV limit. A natural Landau 
pole free form of X{t) is 



In" 



ci + c e 



(19) 
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with a new parameter c > 0. We will see below that cq and ci are fixed by 
the required IR behavior in terms of Ay, a and e^^. c is the crossover parameter 
mentioned in the introduction. 

The relation of the A-parameter in the V scheme, Ay, to that in another scheme 
is fixed by a 1-loop calculation [15]. For instance, the relation to the MS scheme 
is given by [16, 17, 18], Ay/A,^ = exp[(31iV - 10nj)/(66iV - 12n/)]. 

If one wants to incorporate three loops, one would write 

X{t) = In' [ci + e«^* e«i-^^''W] , (20) 

and for e-'^'^^*^*) use the form ln''(c2 + c e*), say. Note that the leading 3-loop effects 
are given by the e"i* term. They are included if we simply replace e"i^^^^*^*) by a 
constant c. 

Eq. (18) is a generalization of Richardson's ansatz [19], who writes PoaR^q) — 
1/ ln[l + g^/A^]. His potential, defined as in (4), has the advantage of simplicity, 
containing only the single parameter Ar. However, by the same token, the corre- 
sponding potential can be fitted to describe only the UV or the IR or the crossover 
region, but not all three or even two of them. For potential models mainly the 
crossover region is relevant, explaining why the Richardson potential works so well 
in that context. 



2.2 The Infrared Regime 

To built the IR behavior given by cq. (1) into our ansatz we have to transform 
eq. (1) to momentum space. A priori the Fourier transform of cxr is ill-defined. 
However, the potential corresponding to av{(i) — q~"' can be defined by analytic 
continuation in n to be 

1 

f (-■) = r(„ + i)co3f • 

as long as n is not an odd integer. So (1) corresponds to the small q expansion 

- + ^ + ■■■ ^ (22) 

and matching this to eq. (18) leads after a small amount of algebra to 

For given N, Uf and Ay, the string tension determines cq. ci is then fixed by the 
IR charge. In other words, Cq and Ci parameterize the leading and subleading IR 
behavior, respectively This was of course already obvious from their definition in 
eqs. (18) and (19). 
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2.3 Summary 



The ansatz defined by eqs. (18), (19) and (23) will be referred to as our "standard 
ansatz" . It is consistent with 2-loop QCD in the UV and in the IR reproduces the 
leading and subleading prediction of the string picture. For given N, Uf and we 
have to fit three parameters in Q;y(g): Ay, cr and the crossover parameter c. The 
leading 3-loop effects can be taken into account by adding one more parameter, Oi 
in eq. (20). At this point we fix the IR charge at = leaving a check of this 
value to future high precision studies. 

3 Fitting and Error Analysis 
3.1 Numerical Evaluation of V^2(r) 

It is crucial to have an efficient way of numerically evaluating V^2 (r) if one wants to 
perform fits. Recall first of all, that our approach is not expected to work for r — 0. 
So we can not simply subtract off the r — value of our ansatz when comparing 
with the Wilson loop potential, which is normalized to vanish for r = 0. by such 
an (implicit) subtraction."^ For simplicity we have therefore added a constant in 
eq. (5), after performing the subtraction. [The subtraction is still useful, since then 
the integral has an integrable singularity at g = 0, and does not have to be defined 
by analytic continuation.] 

Since the integrand of Va2{r) is a periodic function in each qi one can quite 
efficiently evaluate the integral as follows (cf. [20, 21]): First introduce a small 
gluon mass, i.e. replace by l/(g^ + m^). Next approximate the integral by a 
sum; schematically 

(where each component of n runs from —L/2 to L/2 — 1 in integer steps) but only 
after applying the following change of variable to the integral: Qi ^ Qi — esing,, 
where for given m, e is determined by the equations e = 1/ coshw, m — tanhw = m. 
Using the Poisson resummation formula and a contour deformation argument one 
sees that this choice of e greatly accelerates the convergence of the sum towards 
the integral as L is increased. Finally, extrapolate to zero gluon mass. For the 

^This normalization is related to the overall constant in the potential (defined, for instance, 
by the behavior of the potential in the IR). Although this constant is physical — it determines 
the absolute energy of heavy quark bound states, for example — it has not yet been possible to 
calculate it in the Wilson loop approach, basically because it is difficult to separate from self- 
energy contributions. Therefore one defines the Wilson loop potential to automatically vanish at 
r = 0. See [1] for a discussion of issues related to the constant in the potential. 
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accuracy desired here, about 10~^ for small r and 10~^ for large r to be safely 
below the Monte Carlo errors, we can dispense with the last step and simply use 
a fixed small mass such as am = 10"*^. 

For small r we find that L = 20 is sufficient. As the components of r increase the 
integral becomes more and more oscillatory, so L must increase correspondingly. 
If one needs to evaluate ^^(i") for many r- values, as we do, one should make sure 
to evaluate the r-independent part of the integrand only once for a given 27rn/aL, 
instead of doing so again and again for different r. 

For sufficiently large r (depending on the accuracy of the potential data; typ- 
ically when at least one component of r is greater than 16a) one can simply use 
the continuum V{r) given by eq. (4). This integral is divergent in the IR and 
converges slowly in the UV. One can (define and) evaluate it efficiently as follows: 
To improve the convergence in the UV, add and subtract the Richardson potential 
with suitably chosen A^. This potential has a quickly convergent integral repre- 
sentation [19]. Then cancel the singularity of the new integrand by subtracting 
a term with a{q) oc and add it on again in the analytically continued form of 
eq. (21). In the end add a suitable constant to match V{r) to Vo2{v). 

3.2 Fitting Procedure 

As it stands, we would perform fits of the MC data to the ansatz (5) by varying 
the overall constant Vq and the parameters Ay, a and c in av{q)- However, it 
turns out that Ay is strongly correlated with c (and Vq), leading to a large error 
in Ay without a corresponding large error in «y(g). Clearly, the A-parameter is 
a bad parameterization of the physics. It is a much better idea to use ay(g*), 
where q* is some UV scale, as independent fit parameter. For easy comparison 
with the Lepage-Mackenzie method of estimating av{q*) (see sect. 4.2) we choose 
aq* = 3.4018. 

We will therefore fit Vq, ay(g*), a and c. At present all our fits are uncorre- 
lated, i.e. we do not take into account correlations between the potential data at 
different r (we will comment on this in sect. 4.1.2). Instead, we form the usual, 
naive from our ansatz and the Monte Carlo data, which are available to us in the 
form Vmc(i") ± 5Vmc{^)- To find the minimum of we use Powell's method [22], 
which does not require knowledge of any derivatives. 

3.3 Error Analysis 

We will find that the optimal value of c is always quite small, with an error that 
can be of the same order as the average value of c. This is unfortunate, since, we 
recall, for av{q) to be consistent with full 2-loop evolution, any c > 0, but not 
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c = is allowed. The optimal value is somewhat irrelevant, though, because the 
distribution of c is often very non-gaussian, with a long tail towards larger values 
of c. 

We will also find that c can be quite strongly correlated with Vq and ay(g*), 
whose distribution is then also skewed. The reason for this correlation is that, 
in particular for coarse lattices, only the first few points really "see" perturbative 
QCD, and we then have an over-parameterization of the UV and intermediate 
regime in terms of Vq, av{q*) and c. As one would expect, these problems tend to 
go away as the data become more precise and the lattice spacings smaller. 

The correlation between the parameters av{q*), c and a is usually such that 
they tend to move up and down together. The distribution of av{q*) and c can 
usefully, if somewhat simplified, be visualized as a narrow and relatively fiat ridge, 
whose maxiniTim is close to the end of the ridge where both parameters are small. 
In such a situation, error estimates based on the covariance matrix at the optimal 
parameter values are very misleading.^ 

To obtain reliable error estimates in all cases we have decided to use the fol- 
lowing procedure. Once the optimal values of the fit parameters are known, we 
sample the total distribution around these values by varying the fit parameters, 
accumulating all quantities of interest by weighing each sample volume dV with a 
relative probability rfV" exp(— [Confer the maximum likelihood justification 
of the method.] For the fit parameters themselves the number of distinct values 
sampled was only 10 — 20, so to estimate the 1-sigma band and similar statistical 
measures from the corresponding histograms we suitably interpolated the distri- 
bution between these values. For other quantities, like Ay or tq (cf. sect. 5.2), 
the number of different values accumulated is very large and no interpolation is 
necessary. We have checked with various toy examples that this procedure gives 
reliable error estimates. In the long run one should of course perform the error 
analysis by bootstrapping the potential data — an option we do not have at this 
point. 

For the error analysis it would of course be nice if we could replace c by a more 
well-behaved fit parameter for the intermediate regime, namely one that is less 
correlated with av{q*) and Vq and whose distribution is more gaussian. We have 
tried various options — e.g. replacing c by go or qo/^/a of sect. 5, or by Q;y(g*/3) 

^It does not help, by the way, to consider Inc instead of c as fit parameter. Note, first of all, 
that this would correspond to a different assumption about the prior distribution of c (in the 
sense of Bayes' theorem) — now all values of In c, not c, are considered equally likely a priori. 
In this case the distribution of Inc is strongly skewed in the opposite direction, with a long 
tail towards negative values of Inc and a sharp fall-off in the positive direction. In some cases 
the distribution even seems to become unnormalizablc, since the tail to the left does not vanish 
sufficiently fast. We therefore stick to fitting c, using the Inc fits as one measure of systematic 
errors. 
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— but with the present data all have problems of one sort or another in at least 
one case. It is possible though, we think, that some of these problems will not 
exist for more accurate potential data.^ 

The above concerns the statistical errors. Since systematic errors affect different 
quantities to different extents, they will be discussed in the next two sections where 
we present explicit results for various quantities. 

4 Results and Comparison with Other Methods 
4.1 Overview 

We applied our scheme to potential data obtained by Monte Carlo simulations of 
Wilson loops in the following lattice gauge theories: 

• Pure SU(3) at f3 = 6.0, 6.4 [7] on a 32^ lattice and at /? = 6.8 [23, 24] on 
a 48^ X 64 lattice. Potential data exist, basically, at all lattice points that 
are multiples of (1, 0, 0), (1, 1, 0), (1, 1, 1), (2, 1, 0), (2, 1, 1) or (2, 2, 1), and 
have no component larger than 16 for /3 = 6.0 and 6.4, and no component 
larger than 24 for (3 — 6.8. For orientation we jump ahead and remark 
that the lattice spacings of these theories are between roughly 0.10 fm and 
0.038 fm. We should mention that the (3 = 6.8 data from [7], which we 
had first analyzed, turned out to suffer from large finite-size effects — even 
though they are nicely linear at large r — as the comparison with the new 
data shows. We will see later that the new P — 6.S data also show signs 
(much smaller thoTigh) of finite-size and/or other systematic errors. Finally, 
we should point out that the f] = 6.8 data used in [13] are from the same 
Wilson loop data, but used a preliminary extraction [23] of the potential 
from these Wilson loops. New SU(3) data will be analyzed with our method 
in [25]. 

• SU(3) with two flavors of of dynamical staggered fermions of mass am = 0.01 
at /3 = 5.6 [26], and two Wilson fermions (k = 0.1675) at /? = 5.3 [26]. The 
lattice size is 16^ x 32 in both cases, with a spacing of about 0.10 fm and 
0.14 fm, respectively. The potential was calculated at all points in a plane 
with no component larger than 8. There was no sign of string breaking for 
the r-values considered. It seems that the main effect of the non-zero fermion 

^For example, if one uses go, the problem is that the mapping from av{q*) and qo to Ay 
and c becomes nearly degenerate when c is very small. The fact that in some cases the optimal 
value of c is very small, might be an artifact of present data that will disappear with better data 
(more precisely, once the effective potentials have been extrapolated to large euclidean times, 
cf. sect. 4). 



11 



masses is to renormalize the fit parameters. [We will see however that the 
efi^ect of unquenching can not be absorbed in a "/3-shift" with respect to 
quenched QCD.] 

• Pure SU(2) at /9 = 2.85 [6] on a 48^ x 56 lattice. Potential data''' were available 
at all points with no component larger than 3 and at all even on-axis points 
r/a = 4, 6, 8, . . . , 24. The lattice spacing is about 0.028 fm. 

Recall that P is related to the bare lattice couplings go, respectively cto, by /3 = 
2N/g^ = N/27iao. 

4.1.1 Results from the Standard Ansatz 

Using our standard ansatz of sect. 2 we could fit all data down to one lattice 
spacing, with the uncorrelated ~ -^df- In sect. 4.1.2 we will discuss in detail 
what possible effect the neglect of correlations has on our results. To convince 
the reader that our fits really are much better than previous ones. let us here just 
note (more details will be given as we proceed) that in contrast to our ansatz, 
it is completely impossible, at least on fine lattices, to describe the data within 
errors down to one, or even several lattice spacings with the modified Coulomb + 
linear ansatze discussed in sect. 1. This is particularly obvious from plots of the 
difference between fit and data; we will present such plots below. We also note 
that our fit parameters are essentially stable with respect to variations of the fit 
range, which can not be said about Coulomb + linear fits on fine lattices. 

Our results for the optimal parameters and of the fits are given in table 1. As 
already mentioned in sect. 3.3, the distribution of c is quite non-gaussian, and to 
obtain reliable errors we use the procedure described there. The results are shown 
in table 2. As it turns out that the 1-sigma band (determined by cutting off 16% 
of the distribution on each side) is typically rather symmetric around the average 
of a random variable, even if the distribution itself is quite skewed, we saw no need 
to quote asymmetric errors. [In some cases we slightly shifted the mean, to put it 
in the middle of the 1-sigma band.] Note however that the optimal parameters of 
table 1 usually lie very asymmetrically within the 1-sigma band quoted in table 2, 
for av{q*), Ay and c sometimes even outside it! 

Note that in the above tables we also included data sets denoted by /? = 6.8' 
and 2.85'. This was done to illustrate the effect that uncertainties of the Monte 

''As in the references above, for the SU(2) case and those with dynamical fermions we actually 
used data for the "effective potential" Vrir) obtained from Wilson loops at euclidean times T 
and r+ 1, with T = 3 for the Wilson fermion case, and T = 4 for the SU(2) and staggered fermion 
cases. The actual potentials, obtained by a large T extrapolation, will differ somewhat from 
these effective potentials. 
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Table 1: Optimal fit parameters and uncorrelated for our standard ansatz. W, 
S indicates Wilson, respectively, staggered fermions. The (3 = 6.8' and 2.85' data 
sets are explained in the main text. Note that Ay is not a fit parameter in our 
procedure; it is quoted for reference purposes. 
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Table 2: Average fit parameters and 1-sigma band for our standard ansatz, ob- 
tained with the error analysis procedure of sect. 3.3. 



Carlo potential data have on our fit parameters. For the /3 = 6.8' SU(3) data 
set we omitted data points at large distances, namely all points with at least one 
component larger than 16a. This gives us an idea if there are finite-size (or other 
large r) artifacts. These seem to exist, since the string tension from these two 
fits are not consistent. These effects are probably not quite as large as table 2 
indicates, since by incorporating 3-loop effects we will later find better fits, where 
the difference is less significant. But the difference is still larger than in all other 
cases, where the string tension is completely stable when varying the fit range over 
comparable scales (in physical units). 

For the (3 = 2.85' SU(2) set, we increased the errors of the potential at the first 
two lattice points from 1 • 10~"^ to 2 • 10~"^. From the tables in ref. [6] the extrapola- 
tion of the effective potential to large times might easily change the values by this 
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amount, and, in contrast to changes at large r, that presumably can be absorbed 
into a change of the string tension without otherwise affecting the fit, it is not 
obvious that an analogous statement will be true for short distances, where errors 
are much smaller. In this case, increasing the errors at small r significantly reduces 
X^, but docs not lead to much different values for the fit parameters. We regard 
the two "primed" data sets as somewhat more reliable than the corresponding 
"unprimed" ones. 

One might wonder why the fractional error of av{q*), although quite small, is 
much larger than that of the lattice potential at r = a, which in the worst case is a 
few times 10^^. The main reason for this is that ay{q*) is strongly correlated with 
the overall constant Vq of the potential; an increase, say, of Vq can be partially 
compensated by a decrease of av{q*) (and smaller changes of c and a). If the 
potential data at small r are correlated, fits that take this into account should give 
a more precise av{q*) (see below). 

Note the relatively large statistical error of Ay in table 2, to which, as we will 
see below, we have to add a systematic error that in some cases is even larger. 
This large error is not just due to the fact that only the logarithm of Ay enters in 
av{q)- Instead, the main contribution to its error comes from the correlation with 
c. Since both parameters are, roughly speaking, intermediate range parameters, 
one should not be surprised that their correlation is such that it leads to a much 
smaller error in the running coupling than expected from naive error propagation. 
This was the reason we chose av{q*), not Ay, as the fit parameter for the UV 
region. 

Prom table 2 we see that the quenched /3 = 6.0 and the unquenched P — 5.6 
SU(3) theory have almost the same string tension. In previous studies it has often 
been found empirically that the effect of unquenching can be absorbed into such 
a change of /5. Of course, strictly this can not be true, since it would contradict 
QCD at short distances. And indeed, we find in table 2 that these two theories 
have very different UV couplings av{q*)- 

4.1.2 Correlated versus Uncorrelated Fits, Systematic Errors 

We should comment on our use of naive, uncorrelated fits (cf. sect. 3), that ignore 
any correlations the potential data might have between different r-values. Of 
course, one should eventually take these into account, preferably by bootstrapping. 
We do not have this option at present, but would like to dispel any worries one 
might have, that we are, say, underestimating the errors because of this. 

To investigate this question we have performed fits of all potential data to 
modified Coulomb + linear ansatze (cf. sect. 1) using the uncorrelated x^, and 
compared them to previous results from correlated fits with the same ansatz and 
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fit range [6, 7, 26, 23, 27]. We estimate the errors from the covariance matrix, 
which for these fits should be reasonably accurate, since in contrast to the crossover 
parameter c in our ansatz, for example, none of the parameters here is expected 
to be strongly non-gaussian (a posteriori this is also justified by the excellent 
agreement with errors from bootstrapping for weakly correlated data). We find: 

(i) The uncorrelated is smaller than the correlated one; slightly if the data 
are weakly correlated, much smaller if the data have large correlations. 

(ii) The central values of the fit parameters are largely unaffected, even if the 
correlations are strong. 

(iii) The errors from uncorrelated fits tend to be too large, significantly so if the 
data are strongly correlated. 

Basically, these features are well known [23, 27, 28]. In fact, thinking about how 
is affected by (positive) correlations, and how this in turn affects the "Ax^ = 1" 
method of estimating errors, the above features are more or less obvious. The first 
point, that our fits will underestimate if the data are correlated, is of no big 
concern. Our aim in this paper is (a) to show that our ansatz for the potential is 
much better than Coulomb + linear ansatze, and (b) give an estimate of the errors 
of the running coupling extracted from our fits. If point (a) is not clear already on 
purely theoretical grounds — assuming QCD is correct — then one only needs to 
look at some pictures of the fractional difference between our fit and the data (see 
below) and compare them with corresponding pictures from Coulomb + linear fits 
(see e.g. [7]). 

For the cases at hand, the only potential data that have strong correlations 
are those for pure SU(2). In this case, our errors for the parameters of ay(g) are 
likely to be over- not underestimates, according to (iii). [In cases where the first 
few r-values had to be omitted from the Coulomb + linear fits, we can not exclude 
that the potential data have non-negligible correlations among these points. If so, 
correlated fits with our ansatz would allow us to determine av{q*) more precisely.] 

We add some remarks from our experience of comparing (uncorrelated) fits 
to our ansatz and to modified Coulomb + linear ansatze. We found that it is 
somewhat of an art to extract parameters from Coulomb + linear fits, in particular 
on finer lattices. Even with the lattice artifacts taken into account as described in 
sect. 1, and the Coulomb charge e replaced by the mock running charge e — //r, 
one has to omit many small r points from the fit, and it is often not so clear if one 
ever reaches a plateau at which the parameters stabilize. For the /3 = 6.8 SU(3) 
data, for instance, neither e, /, or the coefficient of the lattice correction seems to 
stabilize; / even becomes negative at some point. One clearly sees that this type 
of ansatz does not incorporate the correct physics. 
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Table 3: Optimal fit parameters and for 7^ 1 ansatze. 



Out fits, in contrast, are stable with respect to varying the fit range. ^ Further- 
more, we can at least argue — and we will in sect. 4.2 — that our UV parameter 
av{q*) is most accurately determined at the shortest distances, but no such argu- 
ment is possible for the unphysical "running charge" e — j jr. 

We now turn to the systematic errors of a and c in our ansatz. There 

are two potential sources of such errors: (i) our ansatz for the continuum running 
coupling, and (ii) remaining lattice artifacts not taken into account in our ansatz 
for the lattice potential in terms of the continuum coupling, eq. (5). Of course, 
in general we can not rigorously disentangle and separately investigate these two 
effects, since our fit necessarily involves both. [We can disentangle them at the 
1-loop level, as we will discuss in sect. 4.3.] It is pretty clear, however, that lattice 
effects will have a very small direct influence on a and probably also on c. They 
might have an effect on ay{(f)^ and through the correlations of the fit parameters 
this would then indirectly infiuence c and to a smaller extent also a. 

The question whether Q;y(g*) has errors due to lattice artifacts is sufficiently 
important to be discussed separately in the next subsection. Here we will discuss 
the errors of the fit parameters due to our ansatz for the coupling. Our fits in 
themselves give no indication of significant systematic errors in this ansatz. 

The only real check of our ansatz for Q;y(g), then, is to try to find other ansatze 
that give similarly good or better fits. This is not so easy. Prom sect. 2 we know 
how to incorporate higher loops into ay(g). However, each loop adds at least two 
new parameters (one corresponding to the coefficient in the /3-function, the other 
is a non-perturbative parameter arising from our method of avoiding the Landau 
pole), and this often leads to stability problems with the fits, since the limited 

^There is a slight but noticeable dependence on the fit range for the /3 = 6.8 SU(3) data, where 
(of. the previous remarks on the /3 = 6.8' data set, and below) it is just as likely to be due to the 
data, rather than our ansatz. 
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Table 4: Average fit parameters and 1-sigma band for oi 7^ 1 ansatze. 



part of the intermediate and UV region available in practice has effectively been 
vastly ovcrparamctcrizcd. Another way of saying this, is that with present data 
non-perturbative and all-order effects apparently become important before we can 
clearly discern 3-loop effects. 

We have tried various ansatze that fully or partially take into account 3- and 4- 
loop effects. Even if they give stable results, most give significantly better fits only 
for one of the theories we are considering. This might simply reflect a (not even 
large) fluke of the data in question. The overall best method we found has one new 
parameter ai, that takes into account the leading 3-loop effects. It corresponds 
to the choice X{t) — ln''(ci -|- c e^^*) in eq. (18). Even these fits are not always 
stable and so in tables 3 and 4 we show results for suitably chosen fixed ai. The 
values of are significantly better for the last four theories shown. None of the 
other "sporadically" better fits we found seems to lead to a much larger variation 
in fit parameters (or other quantities obtained from av{q), cf. sect. 5) than that 
observed when comparing tables 1 and 2 with tables 3 and 4. 

Comparing tables 2 and 4 we see that the systematic errors of av{q*) and a are 
never much larger, if larger at all, than their statistical errors. It is significantly 
larger in some cases for Ay- Since c is really a different parameter for different ai, 
it is not surprising that it can also change significantly. More relevant, however, 
is the question of how much av{q) from different fits differ in the intermediate 
region. We find, not surprisingly, that in this region av{q) is typically even better 
determined than at its "edges" . 

For illustration, one of our fits is shown in figure 1. One sees that every little 
bump and wiggle of the data is reproduced by our fit. Since the errors are too 
small to be resolved in graphs of the potential itself, we show in figures 2 and 3 
examples of the fractional difference between our fit and the Monte Carlo data. 

Though these plots do not in any obvious way indicate that there are lattice 
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Figure 1: Comparison of Monte Carlo data (circles with error bars) for the lattice 
potential in pure SU(3) gauge theory at P — 6.0 with ^^(r) from our standard fit 
(lines), calculated on the lattice points and connected with straight lines to guide 
the eye. Shown is an overview (a) and the small r region (b). 



18 



0.004 



V„2 



-0.004 



1 1 




1 




1 






1 1 




1 




1 





3 4 

r/a 



Figure 2: Fractional difference between our ai = f .5 fit and the Monte Carlo results 
for the /9 = 2.85' pure SU(2) potential in the small r region. 
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Figure 3: Same as figure 2, for our ai = 1.5 fit of the potential ol (3 — 5.6 SU(3) 
with two staggered fermions. 
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artifacts not reproduced by our ansatz,^ one could take the fact that the small 
r region contributes a relatively larger fraction to the total as a hint of such 
a systematic error. However, this is at present just as likely to come from slight 
underestimates of the error or shifts in the central values of the potential data 
(compared to the true values) than from systematic errors in our ansatz. We recall 
that obtaining the true potential involves an extrapolation to large euclidean times. 
Uncertainties in the data arise, because either this extrapolation is not done and 
one simply quotes an effective potential from a fixed time, or, if the extrapolation 
is done, it is hard to independently assess its reliability for small r (where the 
quoted errors are very small), since there was no theoretical expectation of what 
exactly the answer should be. This is underscored by the fact that slightly different 
procedures of extracting the potential from Wilson loops tend to lead to correlated 
changes of the potential. This can have a large effect on the of our fit, and lead 
to significant changes of the fit parameters. 

From tables 1 to 4 it seems the theories under consideration fall into two classes: 
Pure SU(3), and the rest. For the former better fits arc obtained for ai < 1 instead 
of ai = 1, whereas for the latter ai > 1 is preferable. Presumably related is the fact 
that in the pure SU(3) cases c is consistently larger than for the other theories. 

The optimistic view would be that this is due to the size of the third /3-function 
coefficient. It is not known in the V scheme, but if the MS scheme is any guide, 
this can certainly not be sole explanation, since in this scheme [29] ai ~ 0.600 is 
the same for pure SU(2) and SU(3), and ai ~ 0.625 for SU(3) with n/ = 2. 

In certain respects the pure SU(3) data seem not to behave as "nice" as the 
others. For example, the of the fits are not as good, and various quantities 
are less stable with respect to changing the form of our ansatz for av{q) or when 
leaving out small r data points (cf. also sect. 5). The latter is true even though 
for pure SU(3) the decrease of when allowing ai 7^ 1 is relatively smaller than 
for the other theories. One reason that might at least partially explain the relative 
stability and other "nice" features of the SU(2) data, as compared with the P — 6.S 
data which naively would be expected to behave quite similar, is that at large r the 
SU(2) data consist only of on-axis points. Another reason could be that the SU(2) 
data have significant correlations in r, at least for large distances [6]. However, 

^In the process of searching for possible systematic errors due to lattice artifacts, we have also 
checked what happens if we use av{q) instead of ay(g) in eq. (5). It turns out that despite the 
slow running of the coupling this does have a significant effect: It is not possible to get good fits 
at small r; the points r/a = \/2, V3, 2 are impossible to fit simultaneously. It is of course not too 
surprising that one should use aviq) in (5), since the integrand of Va2(r) should be periodic in 
all qi. 

^°We know for a fact that this is true for the existing (3 = 6.8 pure SU(3) data, where Gimnar 
Bali generously provided us with potentials extracted by slightly different procedures from the 
Wilson loop data. 
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this does not distinguish the nf — 2 theories from the other pure SU(3) cases. 

One fact that does differentiate the two classes is that only for the pure SU(3) 
data was the potential extrapolated to large euclidean times; in the other cases 
only effective potentials were available. Another point to note is that for the pure 
SU(3) theories potential data exist for much larger distances. We performed a fit 
of the /3 — 6.0 data where we only included points that also exist for unquenched 
SU(3) at P — 5.6. The value of c decreases significantly, and, more generally, the 
SU(3) data behave more similar to those of the second class also in other (though 
not all) respects. 

Although we do not understand why c should decrease, the facts mentioned in 
the last paragraph would seem to provide the best guess, at present, of why the 
theories under consideration should appear to fall in two classes. The "nicer" 
behavior of the second class might disappear once extrapolated to large euclidean 
times. In fact, it was observed in [6] that the difference between the (on-axis) 
effective SU(2) potential and the extrapolated one is not smooth in r. If this is 
the problem, one should investigate if it is possible to perform the extrapolation 
in a smoother way. 

Note that in itself the small value of c in the second class of theories is somewhat 
questionable, since c = is not consistent with full 2-loop evolution. On the one 
hand, it is only In c that enters into the UV expansion, so a small value of c does 
not necessarily violate the principle of naturalness. However, note that the value 
of c is even smaller for the ai > 1 fits than for the ai — 1 ones. This makes a certain 
mathematical sense, but one might question it on physical grounds: Recall that 
the of these fits is significantly better than those with ai = 1. If we take this 
to mean that we are really seeing 3-loop effects, then it is surprising that at the 
same time the 2-loop effects are suppressed due to the small value of c. However, 
since a small c always suppresses both 2-loop and (the leading) 3-loop effects, this 
argument is not conclusive. Only new data, with careful consideration of finite-size 
and extrapolation effects, will help to clarify these questions. 

4.1.3 Estimate of the String Tension 

For future reference we would like combine our results and quote estimates of the 
string tension for the various theories we studied. To obtain the best values we 
should not just take into account the results in tables 2 and 4, but also the following: 
We are performing a global fit and do not want any potential uncertainties at 
small r to distort the fit at large r. Therefore we fix av{q*) at its average value 
obtained from a fit with all data points, and leave out successive short-distance 

^^To be sure, note that the "smoother" behavior of the second class is not due to smaller errors. 
Despite the extrapolation, it is the pure SU(3) data that have much smaller errors. 
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Table 5: Comparison of our estimate of the string tension in lattice units with 
previous ones from (modified) Coulomb + linear fits. 

points until a stabilizes within errors. In most cases this leads to no significant 
change. Combining all information in a semi-objective manner — being more 
systematic about this seems pointless in view of the uncertainties discussed above 
— we obtain the estimates in table 5, where they are compared with previous 
results from Coulomb + linear fits. The final results for the pure SU(2) and SU(3) 
{13 < 6.4) data are quoted in [30], respectively [31]; cf. [23] for /3 = 6.8. The n/ = 2 
results are from [26]. 

As expected, our results are generally in good agreement with previous esti- 
mates, though often more precise. An exception is j3 = 6A pure SU(3), where our 
result is significantly higher than the latest Coulomb + linear estimate (though it 
agrees much better with a previously published estimate [7] and our own estimate 
from uncorrelated Coulomb -|- linear fits). 

Since our method takes into account the corrections to the asymptotic linear 
behavior more accurately, one would expect that it gives better estimates of the 
string tension than previously possible. That our errors are not always smaller 
than previously quoted ones, is due to the fact that we were rather careful about 
systematic errors. In particular, for /? = 6.8 the quoted error is almost completely of 
systematic origin. If one were similarly careful about systematic errors for Coulomb 
-I- hnear fits, their errors would be as least as large, often much larger. 

4.2 Lattice Artifacts, av{q*), and the Lepage-Mackenzie 
Method 

Our precise estimate of av{q*) depends on including all points down to r = a in 
our fits. Leaving out r = a does in some cases lead to significant shifts in the 
central value. This is illustrated in table 6, where we show results from fits where 
data points at small r have been ignored. 

If one had very precise data on a relatively fine lattice, leaving out r = a should 
not prevent an accurate determination of av{q*)- This seems to be the case for 
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Table 6: Optimal fit parameters for some Oi 7^ 1 fits with small r points left out. 
rmin indicates the smallest point included in the fit. Although not shown, the errors 
of av{q*) and especially c rise rapidly as rmin is increased, except for the (3 — 6.8 
data. 



the SU(2) data, as table 6 indicates. However, if the data are not that precise, 
leaving out r = a is expected to lead to much larger errors, for the following reason: 
We will see in sect. 5 that the momentum scale q* corresponds to a distance scale 
of roughly r* ~ 1.53/g* ~ 0.45a. So, ignoring r~a will make it much harder to 
obtain precise information at scale r*. 

One therefore should not be surprised that the determination of av{q*) rapidly 
deteriorates as one leaves out data points on the coarser lattices. A bit surprising 
is the fact that for the (3 = 6.8 SU(3) case the value of ay{q*) from a full fit seems 
incompatible (on the 3 — 4 sigma level) with that from fits without r = a, which 
still have very small errors and a much better (for the (3 = 6.8 results shown in 
table 6 the error of av{q*) is always less than 4 • 10~^). 

Since it is not obvious that our ansatz (5) really incorporates lattice artifacts 
with sufficient accuracy at one lattice spacing, the question arises whether our 
method indeed breaks down at r = a in this case (then the value of av{q*) from 
fits ignoring r = a would be more accurate), or if instead the value of CKv{q*) from 
the full fit is more reliable, and the change when leaving out the first point is just 
the kind of "edge effect" that one might expect if one has data with very small 
(perhaps slightly underestimated) errors, and considers a quantity that depends 
mainly on the "boundary" of the data set. 

Although this question can not really be answered with certainty at present, 
we would like to argue that the second possibility is more likely, i.e. the value of 
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av{q*) from the full fit is quite reliable also for SU(3) ai /3 — 6.8. We already 
mentioned that there are some uncertainties in the [3 = 6.8 data at large r, so the 
same might be true at small r (where previously no independent check of the data 
was possible), in particular in view of the fact that the piuc SU(2) data do not 
seem to have this problem. The lattice spacing in this case is comparable to that 
of the SU(3) theory, even smaller, and the errors at short distances are about as 
good. [At r = a, V2a the potential in the SU(2) case has relative errors of 4 ■ 10~^, 
3 ■ 10^^, respectively, whereas for SU(3) the corresponding numbers are 2 • 10~^, 
4- 10-^] 

These arguments would be more convincing if we had a completely independent 
estimate of av{q*) that we could compare with. Fortunately, such an estimate 
exists. Namely, Lepage and Mackenzie [10] have argued that the running coupling 
av{q) at some UV scale can be determined by comparing the non-perturbative 
expectation value of the plaquette operator, Wu = ;^(Tr f/piaq), as measured in 
the Monte Carlo simulation, with its expansion in terms of av{q)- From the results 
of [16, 17, 18, 11, 27] this expansion can be written as 



In 1^11 = TrCpav{q) 



, JllAT, /6.7117n2 1 n / X2N 

+ 0{ccy{qf) , (25) 



where 5j = — 0.060 for staggered and 5j = — 0.105 for Wilson (r = l) fermions. 

The next question is what q to choose in eq. (25). Were the expansion known to 
large orders, one could extract OLy{(i) at basically any q. For a given truncation of 
the series, however, there is an "optimal scale" a^* at which to extract ay{(i). Intu- 
itively (f is the "natural scale" of the "physical process" in question. Let us write 
the leading non-trivial contribution to a lattice quantity as ayi^q*) /if '^^^ -^(O')- 
[Depending on the quantity in question this might be the 1-loop or the tree level 
contribution.] To this order the natural definition of is^^ 



«y(g*) r''^ d^q F{q) = f''^ d'q av{q) F{q) . (26) 

J— vr/a J—Tr/a 



Using the 1-loop expansion of av{q) around q = q* leads to an explicit formula 
for q*, 



d'qFjq) In g^ 



= ./a . ■ (27) 



According to our philosophy it would seem more natural to use av{q) in eq. (26) to define 
a quantity q*. If this is done, one would also have to use q in eq. (25), for example. However, in 
that equation one wants to think of ay (q) as a continuum object, so the use of q might not appear 
natural. We will not pursue this question here, since in the context of the scale determination 
the use of g or g will not make much difference in the final results. 
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Table 7: The expectation value of the plaquette, the coupling at the UV scale 
aq* = 3.4018 determined by the Lepage- Mackenzie prescription, and the value of 
the A-parameter obtained by assuming ay(g) to be given by the leading part of 
2-loop evolution. 

The integrands F{q) appropriate for planar Wilson loops can be found, for example, 
in [11]. For the plaquette Wn or Inl^n one gets the right q* if one simply uses 
F{q) = const. Using the methods of sect. 3.1 we obtain aq* — 3.4018. 

Expressed in terms of the physical coupling av(g*) the expansion of Inl^n is 
expected to be well-behaved, i.e., at least the first few higher order coefficients in 
the square bracket of eq. (25) are expected to roughly follow the pattern established 
by the first two, which in this case both are 0{1). Using the measured Wu to 
determine av{q*) from (25) with the higher order coefficients set to 0, we will 
therefore quote it with nominal error av{q*)^, which would be the actual error 
(for small av{q*)) if the next coefficient had magnitude 1.^^ The value of av{q*) 
obtained in this way will be denoted by av{q*, Wu). These numbers are collected 
in table 7 for the theories under consideration. 

Note that other quantities besides Wu can be used in the Lepage-Mackenzie 
prescription (see figure 4 below). Each quantity comes with its own q*. What 
distinguishes Wu is that among all simple quantities whose av{q) expansion is 
known to second order it has the largest q*. It should therefore give the most 
reliable results for the coupling in the UV. 

By assuming 2-loop evolution to be accurate at q*, one can obtain an estimate 
of the A-parameter from av{q*)- For future reference we quote in table 7 the 
value of Ay obtained from av{q*,Wu) by assuming that av{q) is given by the 
leading part of 2-loop evolution, Poav{q) — l/{t + blnt). This value is denoted by 

Comparing with tables 1 to 4 we see that our and the Lepage- Mackenzie es- 

^•^ Actually, for gauge group SU(2) the second coefficient is smaller by a factor 2/3 (for n/ = 0) 
when compared to SU(3), and we indeed find that various errors in the SU(2) case are smaller 
by roughly a factor (2/3)^. 
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9.0 
18.0 


0.75614 
0.88450 


0.0731 (4) 
0.03040 (3) 


0.0733 (4) 
0.03049 (11) 


2.4/5 
3.1/5 


0.0747 (12) 
0.0311 (4) 


0.4/4 
0.8/4 



Table 8: Comparison of av{q*) for pure SU(3) from our fits (fourth and sixth 
columns) with results from the Lepage- Mackenzie prescription. The latter is quoted 
with error av{q*Y- 



timate of av{q*) are consistent within their errors, our values being somewhat 
smaller (in particular for theories on coarser lattices). It would be nice to compare 
our and the Lepage-Mackenzie estimate in a situation where the errors are much 
smaller. This is indeed possible. In [32] Wilson loops up to size 8x8 were gener- 
ated for pure SU(3) on very fine 16^ lattices. We used their P — 9.0 and 18.0 data 
data to extract the potential (using times T = 7 and 8) at the first seven on-axis 
points, and then fitted them to our standard ansatz. Obviously a can not be fitted 
from such data, so we simply set it to some reasonably value estimated by scaling, 
the precise value being irrelevant. For simplicity we also fixed c = 0.001 (again, the 
precise value of c has hardly any effect). We performed fits using all r > a data, 
and also fits with the first point left out. In table 8 we show our results for av{q*) 
and the corresponding results obtained by the Lepage-Mackenzie method, which is 
very accurate for these fine lattices. Note the perfect agreement for the r>a fits. 
When leaving out the first point, the central value of the fitted av{q*) becomes 
significantly worse — despite the much better — though with its larger errors 
it is still (barely) consistent with the Lepage-Mackenzie estimate. 

A significant change of av{q*) combined with a much better value for when 
leaving out the first data point — this was exactly the kind of behavior that 
raised some doubts about our estimate of av{q*) for the SU(3) theory at P = 6.8. 
But table 8 suggests that our estimate of av{q*) from fits with all data points is 
reliable. Before we can draw this conclusion with certainty, it will, however, be 
necessary to apply our method to more, in particular more accurate data. In the 
next subsection we perform one more important check of our method; it allows us to 
separately test our ansatz for the lattice potential, independent of any assumption 
about the form of the running coupling beyond perturbation theory. 

To conclude this subsection, we mention that we investigated another aspect 

^''Note that at fixed r/a fractional lattice artifacts do not vanish as the lattice gets finer. 
So there is no obvious reason why our method should work well on fine lattices — meaning 
not just that wc get a good fit of the potential data, but also, that this involves no significant 
compensating error between our handling of the lattice artifacts and the ay {q) extracted from 
the fit — but not very well on coarser ones (it should work somewhat better on fine lattices, since 
the coupling is smaller). 
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of the Lepage- Mackenzie prescription, namely to what extent av{q*) agrees with 
{27:)~^ J^^Jj^d^q av{q)- Using av{q*) and av{q) from our fit, we find that the 
integral is always slightly larger, but the difference is well within the expected 
av{q*)^ error (the difference is never more than half this error). 



4.3 The 1-Loop Expansion of the Potential 



We woTild now like to compare the 1-loop expansion of our ansatz for the lat- 
tice potential with the 1-loop results for the on-axis lattice potential that can be 
obtained from the more general results for planar Wilson loops in finite volume 
presented in [11]. 

Expanding the coupling of a generic scheme around q — /i to one loop gives 
a{q) — — Q;(//)^/3olng^///^. Inserting this into eq. (5)^^ we obtain 



Va2{r) = a(/i)K^[l](r) + a{fifPoVa4^fiVq']{r) + 0{a' 



(28) 



where we use the notation [a] to denote the "potential" in eq. (5) as a functional 
of the "running coupling" a{q). For any given a{q) these functions can easily be 
evaluated numerically as described in sect. 3.1. Remember that Va2[l](r) is the 
lattice Coulomb potential. The expansion in (28) can easily be extended to higher 
orders. 

On the other hand, the 1-loop expansion of the on-axis potential given in [11] in 
terms of the bare coupling can be reexprcssed in terms of a continuum coupling. 
For the latter we choose the V scheme. Restricting ourselves to pure SU(iV) from 
now on, we have (cf. the references given above eq. (25)) 



1 - avil^)(Poln(—Y + 2.409830 TV - 



+ 0{a') , (29) 



which allows us to rewrite the on-axis potential of [11] (in the infinite volume limit) 



as 



V^P{r) = av{fi)Va2[l]{r)+av{fiy 



(lolgs)'^"^^^^^'^ + 

Here, in terms of the notation of Heller and Karsch [11], 



X{R) = (47r)2 hm X{R,T)-X{R,T-1) + - 



W2{R,TY-W2{R,T-lf 



(30) 



(31) 



^^Where one should ignore the constant Vq in the following. In the end we will consider 
potential differences anyhow. 
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18.0 
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0.02343 
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SU(2) 
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0.234 


0.0941(2) 


0.0951 


0.0932 


0.0895 



Table 9: Results for the difference of the lattice potential in pure SU(3) and SU(2) 
at the first two on-axis points from Monte Carlo simulations compared with per- 
turbative expansions in ckv{(1a)j where ag^ = 1.028. The last column contains the 
1-loop result of [11]. 



which is independent of A^. 

Let us now consider the difference AVa^ = aVa^ {{2a, 0,0)) — aVa2{{a,0,0)) of 
the lattice potential at the first two on-axis points. For this quantity the optimal 
choice of fj, according to the Lepage-Mackenzie prescription is /j, — 1.028/ a, which 
we will denote by q\. Using this value, our ansatz gives 

AVa2 = Cf av{qX) [0.54255 - 0.08186 /5o av{qX)] + 0{a^) . (32) 

From ref. [11] and unpublished results of Heller [27] on 32^ lattices we estimate 
X{2a) - X{a) = 0.331(1), so that 

AV}P = Cf av{ql) [o.54255 - 0.234(7) /3o av{ql)\ + C(a^) . (33) 

In tabic 9 we compare Monte Carlo results for AVa2 with the tree level and 1-loop 
estimate from our ansatz and from ref. [11]. The errors of the Monte Carlo results 
were obtained by quadratically adding those of the individual potential values, 
which is almost certainly an overestimate because of correlations. To keep all 
estimates strictly perturbative, we did not calculate av{q\) from our fitted ay (5), 
but instead used the Lepage- Mackenzie value of av{q^) from the plaquette and 
2-loop evolution. 

Note that the 1-loop estimate from our ansatz is consistently and significantly 
better than the result of [11]. This is somewhat odd; why should the 1-loop result 
from our ansatz be better than what is presumably the exact answer? Probably we 
are just lucky with our ansatz. The difference between the Monte Carlo results and 



^^If we use our fitted (xv{q) to calculate the value of ay(g^) to be inserted in eq. (32), the 
agreement between our 1-loop result and the Monte Carlo data becomes even better. For example, 
in the SU(2) case we then obtain AT42 = 0.0941. 
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those of ref. [11] is consistent with a 2-loop contribution of 1.5(2) Cp Po C({Qa)^ to 
AVa2. The coefficient 1.5(2) seems a bit large, perhaps, but it is a nontrivial fact 
that this value is consistent with all data in table 9. We have checked that the 
contribution to AVa2 from zero- modes [33], which was not included in [11], is 
negligible in the cases considered. 

The situation illustrated in table 9 requires further investigation. But we can 
certainly conclude that our ansatz incorporates the lattice artifacts very well. We 
emphasize that our results in table 9 do not involve any assumption about the 
form of av{(l) beyond one loop. It is therefore a check of our ansatz for the lattice 
potential that is completely independent of our ansatz for the continuum running 
coupling av{q)- 

4.4 Comparing av{q) with 2-Loop Approximations 

We would now like to address the question to what extent the running coupling 
obeys 2-loop evolution. The first question, though, is, which form of 2-loop evolu- 
tion we should use. There is Poa{q) — l/{t + blii{t + blnt)), which incorporates the 
leading and sublcading 2-loop contribution. More commonly used are however the 
expressions j3Qa{q) = l/{t + bint), or l3o<y{q) = (1 — blnt/t)/t, which incorporate 
only the leading part of the 2-loop contribution to the running coupling. We will 
refer to these expressions as the "full" , the "leading" , and the "other leading" form 
of 2-loop evolution of the coupling, respectively, and generically denote them as 

We would like to compare these expressions with our fitted av{q), which, 
though not exact, ought to be much closer to the exact answer than any of the 
above expressions. In table 10 we compare pure SU(3) fits for av{q) with the 
various forms of 2-loop evolution. The A-parameter used in the 2-loop expressions 
is Av from our fit. We used the relative scales of these theories, as determined in 
sect. 5 below, to give a comparison at the same physical momenta for the different 
theories. The momenta are quoted in units of the lattice spacing of the f3 = 6.0 
theory. Using either tq — 0.5 fm (cf. sect 5) or y/a ~ 440 MeV, one sees that the 
momenta chosen are roughly 1 GeV, 2 GeV, 10 GeV and 80 GeV. 

Somewhat surprisingly, the leading form of 2-loop evolution is closer to our 
fitted Oviq) than the full one, except for the highest momenta, even though our 
ay{q) was designed to contain all 2-loop contributions. The reason for this is 
that up to moderate UV momenta (around 10 GeV) there is a partial cancellation 
between the subleading 2-loop and non-perturbative contributions. Note that the 
other form of leading 2-loop evolution always fares worse in our comparison than 
either the leading or the full expression. This is not surprising, since compared to 
the full 2-loop formula it has a subleading 2-loop term of the wrong sign. This 
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0.5 


0.493 
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0.261 


1.137 
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Table 10: Comparison of ay{q) from our standard fit for pure SU(3) (using the 
average values from table 2) with the various forms of 2-loop evolution. Ay from 
our fit is used in the 2-loop formulas. 

form of 2-loop evolution should therefore never be used, and we will ignore it from 
now on. 

While at 2 GeV the difference between our av{q) and 2-loop evolution can be 
10 - 20%, at 10 GeV it has already dropped to 1 - 2%. At 80 GeV the difference 
is only a few permille if one uses the the full 2-loop expression. 

Finally, in figure 4 we present a graphical comparison with 2-loop evolution. 
We also show the values of av{q) obtained at various values of q using the Lepage- 
Mackenzie method for various small Wilson loops and Creutz ratios. The rightmost 
circle corresponds to the value extracted from H^n, which agrees best with our fit, 
as expected. 

If one wants to call the agreement between the 2-loop expressions and our 
non-perturbative running coupling good or bad depends on one's expectations. 
The main point is that we have a non-perturbative, numerically accurate running 
coupling, and can quantify how good 2-loop evolution is. 
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av{q) 




1 10 



aq 

Figure 4: Log-log plot of av{q) for (3 = 6.0 pure SU(3) from our standard fit (the 
solid lines delineating the error), its leading 2- loop approximation (dotted line), 
and results obtained by applying the Lepage-Mackenzie method to various Wilson 
loops and Creutz ratios (circles), shown with error av{q)^- 

5 Scaling 

5.1 Introduction 

By the term scaling we circumscribe three related but distinct issues that one 
faces in simulating a quantum field theory on a lattice. First, choosing and then 
calculating a dimensionful quantity from (results of) the simulation. This quantity 
will then be known in lattice units, and one can use it to determine the relative 
scales of different lattice theories. Secondly, one must determine the absolute scale 
of a lattice theory by equating a dimensionful quantity from the simulation with its 
experimentally measured value. Thirdly, one must extrapolate to the continuum 
limit by simulating the theory on several (sufficiently small) lattice spacings; or at 
least check that one is in the scaling region where dimensionless ratios of physical 
quantities are independent of the lattice spacing within some acceptable error. 

The absolute scale determination is of course not really meaningful unless one 
is sufficiently close to the continuum limit — and to the real world, which might 
involve further extrapolations, e.g. in various masses or the number of flavors. 

Note that there is no reason for the relative and the absolute scales to be 
determined from the same quantity. In fact, this is probably not a good idea, 
since it would be rather surprising if the quantity that can be determined most 
precisely from the simulation happens to also be the observable most suitable 
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for precise experimental determination. One should therefore clearly separate the 
absolute from the relative scale determination. Once the latter has been calculated 
for various theories, one just has to determine the absolute scale for one of these 
theories — using perhaps a quite different observable — to know it for all theories. 

In the past it was customary to use the string tension a to set the scale. How- 
ever, since it is an asymptotic (large r) quantity, it has relatively large statistical 
and systematic errors. Furthermore, a is not really well-defined in the presence of 
sea quarks, due to string breaking.^''' Is is therefore also experimentally not par- 
ticularly well determined. The situation is even worse for the A-parametcr, whose 
error is much larger, due to, as we discovered, its correlation with higher order and 
non-perturbative contributions to the running coupling at small and intermediate 
distances. We now describe a better way of setting the scale. 

5.2 The Scales go and ro 

Instead of using a parameter describing the asymptotic long or short distance be- 
havior of the running coupling to set the scale, it is clearly a much better idea to do 
so using the value of the coupling itself at some intermediate distance. Sommer [12] 
proposed to determine a scale ro = ro(cr) using the interquark force, 

Cp OLF{rQ) = Cr , (34) 

for some fixed constant of C(l). [For details about the coupling in the force 
scheme we refer to sect. 6.1.] This is basically equivalent, as we will see below, to 
defining a momentum scale qQ — q^ic,^ via 

Cf av{qo) = Cq , (35) 

for suitably chosen Cg. 

In [12] it was suggested to take Cr = 1.65, since in various phenomenological 
potential models of heavy quarkonium systems this gives ro — 0.5 fm, roughly the 
distance at which these effective potentials are best determined empirically. Since 
the relation between the static potential and these effective potentials is not well 
understood at present, ro does not seem to provide the best way to determine 
the absolute scale. The various extrapolations required to connect Monte Carlo 
data to the real world are much better understood — in particular have been 
explored by simulations — for spin-averaged level splittings in heavy quarkonium 

^^Though in practice the string tension is at present "quite well-defined" , since clear signs of 
string breaking have not been seen yet; instead one observes a large linear regime in the static 
potential (cf. the fermionic theories of sect. 4). That such a regime should exist before the 
potential flattens and keels over is also clear from the success of Coulomb + linear potentials in 
describing heavy quarkonium systems. 
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systems [34, 35], which are of course also accurately known experimentally. When 
available, this approach seems to provide the best method for an absolute scale 
determination (for a review of the controversy over the difference in the absolute 
scale determinations from this method versus those from light quark spectroscopy, 
see [35]). Be that as it may, we will here concentrate on determining relative scales, 
for which a quantity like Tq is ideally suited. 

The scale Tq = ro(1.65) has already been used in the literature, so for compar- 
ison we will calculate it too, using the exact relation between av{q) and a^(r). 
From our point of view, it is however to some extent more natural to define a 
scale qo directly in terms of av{q) using eq. (35). The question is if one can find 
a simple relation between qq and Tq. In practice this is indeed the case. We will 
see in sect. 6.1 that to two loops apir) — av{q) with r = /q — 1.526205/g.-'^^ 
The scales we are interested in are however not in the perturbative region, so this 
equality does not hold very well. Nevertheless, one might suspect that by either 

(i) setting Cq — there is an almost universal relation tq — Crq/qo with some 
Crq 7^ e^"^, or 

(ii) still requiring tq — e^~^ /qo there is an almost universal constant Cg{cr) ^ 
that achieves this. 

Studying the different theories of sect. 4 we find that both of the constants c^q in 
(i) and Cq in (ii) are universal on the permille level for = 1.65. We decided to 
basically follow approach (ii) and define go = go(l-45). With ro = ro(1.65) from 
now on, we find that go = 1-524(4) holds for all theories and all paramcterizations 
of OLy{(i) discussed in sect. 4. [Considering only the standard parameterization of 
Q;v(g) we find go ro = 1.5266(12).] This simple relation between go and ro is useful 
for quick translations between scale determinations in the two schemes. Below we 
will however always calculate ro and go from their original definitions, given by 
eqs. (34) and (35). 

5.3 Results 

Using our standard fit and the exact relation between OLyi^q) and olf^T^ wc obtain 
the results in tabic 11 for go, ro and the dimensionless quantities q^/ y/a and r^^fo. 
These estimates were obtained in the same analysis that lead to the results in 
table 2, including the error estimates from the method described in sect. 3.3. In 
particular, the error of a "composite" quantity like was not calculated by 

^^It is in this sense that a momentum scale q corresponds to a distance scale r ~ 1.53/g, as 
mentioned in sect. 4. Note, however, that one could use other prescriptions to define numerically 
different momentum versus distance scale relations. In the present context the above seems the 
most natural one. 
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Table 11: Average values and errors of the scales go and ro in units of the lat- 
tice spacing and the string tension, obtained from our standard ansatz and error 
analysis. The /9 = 6.8' and 2.85' data sets are explained in sect. 4.1. 
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Table 12: Same as table 11, for ai^l fits. 



naive error propagation from those of ro and cr, which would ignore the correlation 
of these quantities. 

In table 12 we show the corresponding results for our ai 7^ 1 fits (cf. table 4). 
Using the difference to our standard fits as an indicator of systematic errors, we 
see that these are very small for go and ro- We also find that these quantities are 
stable within errors when omitting small r points from the fit. The small amount 
they do move, tends to brings values from different fits closer together. Note that 
the slight discrepancy in go and Tq from the 6.8 versus the 6.8' data sets for ai = l 
disappears for the much better ai = 0.5 fits. 

All this is very good news: Despite the fact that different parameterizations 
incorporating the 2-loop running of the coupling lead to slight systematic errors at 
the "edges" of ^^/(g), that is, in ay(g*) and o", they introduce no sizable systematic 
bias in the intermediate region covered by the potential data. 
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The dimensionless quantities q^sj \fo and r^s/a^ on the other hand, can have 
significant systematic errors, inherited from o. For the SU(2) case these systematic 
errors are still almost negligible, whereas for SU(3) at /9 = 6.8 and the unquenched 
theories the systematic errors are much larger than the statistical ones.^^ 

We combine our results, including those from fits where small r data points have 
been left out, and in table 13 give estimates of ro and that include systematic 
errors. We compare our results with estimates from various other methods: The 
/3 = 6.0 and 6.4 values of ro [31] were obtained by the method of ref. [12] using force 
data from [7, 23]; the value of a used in roy^ is from the same reference. The 
/3 = 6.8 results were obtained with a modified Coulomb + linear fit [23], as were 
the rif — l results [26]. The /? = 2.85 SU(2) value of ro was obtained in a simulation 
using the Schrodinger functional coupling q;sf(?) [36] (cf. sect. 6.2). 

Wc recall (sect. 4.1.2) that our errors for a given fit of the /3 = 2.85 SU(2) data 
might not be accurate due to correlations. However, we saw that our uncorrelated 
fits tend to overestimate errors in such a situation, so all errors in table 13 should 
be on the safe side. 

As another check we have determined ro and roy^ using uncorrelated fits 
to modified Coulomb + linear ansatzc. Despite the uncertainties of the short- 
distance parameters of these fits, in the intermediate range these fits are rather 
well-determined and we obtain results in good agreement with the estimates from 
our ansatz. 

We hope that an analysis of new data at other /3-values for pure SU(3) will 
lead to smaller errors for rov^ [25]. Even as it stands, though, our errors for ro 
and rov^ are much smaller than those obtained previously (except for rov^ in the 
/9 = 6.8 case, where the Coulomb + linear estimate in table 13 does not include all 
systematic errors). Note in particular that for the first time we can see a significant 
difference in ro\/a for theories with and without dynamical fermions. 

As table 13 indicates, for pure SU(3) r-Q^fa scales within errors of about 2% 
for j5 > 6.0 (perhaps even for smaller P which we have not studied so far). Note 
that if we just consider our standard fit results in table 11, ro-v/o" even scales at 
the 0.5% level. Scaling of the string tension is shown in another way in table 14, 
where we used our estimates of a from table 5 and rescaled them to the (3 = 6.0 
lattice spacing with the help of the relative scales obtained from ro in table 13. 

In this table we also show, as check of UV scaling, results for ay (g) at the scale 
5*g = 3.4018/a6.8 for different /?. [A rough check of UV scaling, without error bars, 
was already provided l)y the av{q) entries in table 10.] We used our ai = 0.5 fits 

^^Note, by the way, that the statistical error of rov^ for the last three entries of tables 11 
and 12 is much smaller than for the pure SU(3) cases, even though the opposite is true for the 
errors of the potential data themselves. We had observed and discussed similar features already 
in sect. 4.1 
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Table 13: Comparison of determinations of Tq and r^s/a from our and other meth- 
ods. The references for the latter are given in the main text. 



/5 






6.0 
6.4 
6.8 


0.0491(8) 

0.0496(13) 

0.0519(22) 


0.1114(23) 
0.1164(17) 
0.1142(8) 



Table 14: String tension for pure SU(3) in units of the /3 = 6.0 lattice spacing, and 
Q;y(g) evolved to the common scale ^g g. 

and evolved the /5 = 6.0 and 6.4 results using the appropriate ratios of Tq. The 
error from the scale uncertainty is negligible. Instead, the error is dominated by 
the error of our fit (where one must be careful to take into account the correlation 
of the fit parameters) and by the systematic error of OLy{(f\ which we took to be 
the difference between the central values in tables 2 and 4 (for /3 = 6.8 we averaged 
over the 6.8 and 6.8' results). Wc see that scaling violations, if existent, are just 
slightly larger than the error of 1 — 2%. 

Without presenting explicit results here, we note that scahng holds at the same 
or even shghtly better level of accuracy also in the intermediate regime^° (to some 
extent this can be seen, again, in table 10). We hope to provide more accurate 
results for pure SU(3) in [25], including a more stringent check of the scaling of 
Q;y(g) at all momenta. 

Finally, we have repeatedly stressed the large systematic errors of the A- 
parameter, that are usually ignored. To underscore this point, we compare in 
table 15 three different determinations of this parameter: First the value from our 

^•^One should not be worried about the fact that the crossover parameter c does not scale too 
well (for given a\\ it obviously should not scale when comparing results with different ai). It 
is an effective parameter, representing small subleading effects, and therefore prone to amplify 
even small flukes in the data. 
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Table 15: Comparison of our, the Lepage-Mackenzie, and the estimate of Ay from 
force data using the 2-loop ap^r) [37, 6]. Our value of Ay includes an estimate 
of systematic errors. The last column is our estimate of Ay in units of the string 
tension. 



fit, second that obtained by the Lepage-Mackenzie prescription from the plaquette 
by assuming (leading) 2-loop evolution (cf. table 7), and thirdly that obtained by 
matching force data to the 2-loop apir). We used Ay ~ 1.526 • Ap to translate 
from the force to the V scheme (cf. sect. 6.1). Note that our result is the only 
one taking into account the most important systematic errors, which we estimated 
from tables 1 to 4. Our estimate is systematically lower, but agrees quite well with 
that obtained from the plaquette. As a somewhat "old-fashioned" check of scaling 
we tabulate our estimate of Ay/y^ in the last column. Within large errors scaling 
seems to hold also for this quantity. 

6 Other Schemes for the Running Couphng 

We now briefly discuss two of the more established methods of obtaining the run- 
ning coupling in different schemes. Many other schemes have recently also begun 
to be explored, see e.g. [38, 36, 39]. 

6.1 The Force Scheme 

We already mentioned this scheme at various points in this paper. We now discuss 
it in more detail. Let us first describe the close relation of the coupling in this 
scheme, apir), to Q!y(g). Recall that ap{r) is defined in terms of the continuum 
potential as 

rV(r) = Cpap{r) . (36) 

So just from their definitions we have an explicit exact relation between these two 
couplings. This relation can be rewritten as 

(J 2 f°° dq 

(^F{r) = — H / — singr [qdqav{q) + av{q)\ , (37) 

Cp 71 Jo q 
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where the "IR subtracted" av{q) is defined by 



av{q) 



av{q) 



2a 1 



(38) 



Besides this exact though perhaps not directly iUuminating relation, there is a 
"deeper" relationship between these two couplings. Namely, ^^(r) is in a rather 
literal sense a position space analog of av{q)- To see this, note first of all that, hke 
any continuum coupling, apir) obeys the /3- function equation with the first two 
universal coefficients. In particular, it has an expansion of the form (10), where 
now t = ln(Ai7'r)^^. The relation between Ap and Ay can be obtained by a direct 
calculation of the 1- and 2- loop terms in V (r) , written as a Fourier transform of 
av{q)- We obtain A^? = e'''~^Ay, where 7 = 0.57721566 ... is Euler's constant. 

Next, note that for long distances 



which is to be compared to the IR expansion (22) of av{q)- So we see that one 
can use the same ansatz (18) for aF{r{t)) as for av{q{t)). 

We mentioned in the introduction that the force scheme has been widely used [5, 
6, 7, 37, 40] to obtain an estimate of the A-parameter by matching force data to 
the (leading) 2-loop expression for aplr). In sect. 5 we also discussed how ap{r) 
is used to define the useful intermediate scale Tq by CpaF{rQ) = 1.65. Sommer [12] 
determined tq from force data. Let us note some potential and real problems in 
using the force scheme to determine Ap or ro, and to what extent some of these 
problems can be circumvented. 

• To get the force from potential data involves taking numerical derivatives, 
which generically increases the errors dramatically. Only if the potential data 
at neighboring points are strongly correlated will this not be the case. 

• Lattice artifacts in the force can easily be taken into account at tree level, 
cf. [12]. However, we know that at short distances this is not sufficient to 
describe the lattice artifacts accurately. For the Tq determination this will 
presumably lead to significant errors only, if at all, on rather coarse lattices, 
where rg/a is not much larger than 1, cf. table 13. But for the Ap estimate 
this will induce sizable systematic errors even on fine lattices. 

• Global versus local fits: To determine ro from force data one only needs a 
locally accurate interpolation between two neighboring points [12]. ro can 
therefore be obtained without any assumption about the global from of the 
force. By the same token, however, one misses the opportunity of obtaining 




(39) 
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a more precise value for ro offered by a globally accurate ansatz (or accurate 
at least in the intermediate region) that takes into account information from 
many more data points. For the estimates one often has the worst of 
both worlds, namely, is typically extracted locally from the force using 
the (leading) 2- loop formula for ap{r), which can not be expected to hold 
very accurately. This problem in itself could be mitigated by fitting to a 
global ansatz like (18). However, lattice artifacts would still be a problem. 

Note that none of the above problems occurs in our approach: We have a globally 
accurate^^ ansatz for the potential that takes lattice artifacts into account much 
better than at tree level. What we have not done yet, is to take into account 
correlations of the potential data. If these were strong, so that the force could be 
extracted with great accuracy, the best method to fit the force data would pre- 
sumably be to match them to the corresponding differences of the lattice potential 
as given by our ansatz. This would solve the one problem that our approach has, 
namely that we have to fit the overall constant of the potential, whose correlation 
with av{q*) provides the main contribution to the latter's error. 

However, in most cases the correlations of the potential data do not seem to be 
strong enough to outweigh the other disadvantages of using the force. And, in any 
case, since the potential data are primary even in the presence of correlations, the 
cleanest method of extracting a running coupling would be to append our fitting 
routine to that for the lattice potential and perform fits on bootstrap copies of the 
potential data. This would automatically take care of all correlations, presumably 
in the best way possible. 

6.2 The Schrodinger Functional Coupling 

Liischer et al [41, 30, 31] have introduced and studied in some detail a completely 
different running coupling. It is defined by the response of the partition function 
( "Schrodinger functional" ) of a lattice gauge theory on a finite lattice to a constant 
background field, with the extent of the lattice providing the running scale. This 
coupling is denoted by asF(q'). The great advantages of this method are (i) that the 
coupling is defined so that it can be extracted from the Monte Carlo simulation 
with basically no systematic error and (ii) that it is relatively easy to cover a 
large energy range by varying the size of the lattice. Its disadvantages are (i) 
that it involves "special purpose" simulations due to the gauge field configurations 
required and (ii) that there are apparently large higher order coefficients relating 
this scheme to more common ones, see below (the 2-loop relation to the MS scheme 
will be known soon [42]). 

^^Recall that we checked that leaving out points at small or large distances from the fit does 
not significantly change our estimate of tq. 
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6.3 A Comparison 

We would like to end this section with a quantitative comparison of ap, ay and 
asF at a suitable physical scale. For pure SU(3) at (3 — 6.5 the most accurate value 
oi ap obtained by the UKQCD collaboration from force data [40] turns out to be 
at r/a = 2.5322, where they find q;^ = 0.248(2)(1). On the other hand, translating 
results for q;sf(q') via perturbative matching formulas gives ap = 0.205(7)(9) at 
the same physical scale [31] (we quote ap ~ 0.009 as the second, systematic error 
from the perturbative matching). We do not have results at P — 6.5. But using 
ro/a = 11.23(21) for P — 6.5 (from [31], obtained using data of [40]) and our 
values of vq we obtain Q;i? = 0.231(5), 0.237(5), 0.231(4) at the same physical scale 
for (3 = 6.0, 6.4, 6.8, respectively. In our error estimate we linearly added the 
contributions from our av{q) and those from expressing the above scale r in the 
appropriate lattice units. The systematic error from the ai — 1 versus ai = 0.5 fits 
is negligible. The non-monotonicity of the ccp- values as a function of (3, which we 
previously saw for av{q) (cf. tables 10 and 14), is perhaps due to systematic errors 
of the potential data we used in our fits. We therefore simply quote 0:^ = 0.236(6) 
at the above scale. This is in reasonable agreement with the value from force data, 
but they both differ by about 3a^ with respect to the value from the Schrodinger 
functional. This indicates that perturbative corrections in the translation between 
asF and ay or ap have quite large coefficients. 

7 Conclusions and Outlook 

Our ability to fit lattice potential data to very high precision rests on two ingre- 
dients. First, our expression for the lattice potential in terms of the continuum 
running coupling av{q), eq. (5), and, secondly, our parameterization of av{q) that 
incorporates both the short distance QCD and the long distance flux tube predic- 
tions for the static potential (sect. 2). 

Our ansatz for the lattice potential is obviously correct at tree level. At the 
1-loop level we found (sect. 4.3), rather confusingly, that our ansatz seems to 
work much better than what is presumably to be the exact answer [11]. This 
remains to be understood. Be that as it may, this and other checks, like the 
comparison with the Lepage- Mackenzie estimate of the short-distance coupling 
(sect. 4.2), indicate that we are handling the lattice artifacts extremely well. More 
precisely, the fact that wc can reproduce basically all lattice artifacts within the 
very small errors of the potential data docs not seem to involve any significant 
compensating error between our ansatz for the lattice potential and our estimate 
of the running coupling av{q) at short distances. To be absolutely sure about this 
conclusion, however, we do need new potential data, since at the level of precision 
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we are working many presently available data seem to suffer from slight systematic 
uncertainties of one sort or another (finite-size effects, extrapolation of Wilson loop 
ratios to large euclidean times). 

Previous methods of fitting the data did not incorporate QCD very well — 
to the extent that (slight) systematic errors of the data did not make much of a 
difference; they could not be detected anyhow. Our method is precise enough to 
make it necessary to carefully check the data for finite-size effects and extrapolate 
the effective potentials to large euclidean times. Put in another way, there is a 
large payoff again for precise potential data. Particularly exciting is the prospect 
of potential data from improved actions [43] , for which our method is ideally suited 
(cf. the introduction), since what eventually will prevent us from using our method 
to extract av{q) ever more accurately, are lattice artifacts, and not, for example, 
a 2-loop approximation used in other methods. 

We believe that performing fits of potential data with our method allows for 
the first time a realistic error assessment, in contrast with the hard to estimate 
systematic errors afflicting (modified) Coulomb -|- hnear fits and the method of 
estimating the A-parameter by matching force data to 2-loop formulas. 

We showed that the A-parameter has a relatively large error (cf. sect. 4.1.1) — 
without leading to a large error in av{(l) itself — because A is strongly correlated 
with higher order and non-perturbative contributions at an intermediate range, 
appearing here in the form of the crossover parameter c. This is not surprising, 
if one remembers that the true UV regime is 5 ^ A, i.e. A is really more an 
intermediate range parameter; it must be correlated with higher order effects. We 
therefore regard previous estimates of A as overly optimistic. 

Within our approach the error of A is largely a red herring. We only use A, 
which is clearly a bad way to parameterize the physics, because it is otherwise 
difficult to write down a closed form expression for the running coupling. But the 
physics is contained in av{q), which has much smaller statistical and systematic 
errors. 

The scale tq is now widely used to compare the lattice spacings of different 
theories. For such relative scales, this is certainly much more accurate than using 
the A-parameter or the string tension. Sommer [12] described a careful procedure, 
taking lattice artifacts into account at tree level, to determine Tq locally from 
force data. We determine Vq from our fit of potential data, using the exact relation 
between the force and potential schemes. Our values for tq (sect. 5) are significantly 
more accurate than those from the procedure of ref. [12] — or those obtained 
using a modified Coulomb -|- linear ansatz for the potential — because (a) we use 
potential data, that in the cases considered have much smaller errors than force 
data, (b) we can take lattice artifacts into account more accurately, and (c) we 
have a globally accurate ansatz for the potential, so that information from many 
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points, not just two neighboring ones, can be taken into account. 

Using ro to set the (relative) scale of different lattice theories, enabled us to 
check the scaling of av{q) for pure SU(3), cf. sect. 5.3. If existent, the scaling 
violations are not much larger than the error, which is about 2%. Unfortunately, 
the pure SU(3) data on the finest lattice, corresponding to P = 6.8, show some signs 
of systematic errors, making a more accurate check of scaling impossible, since we 
investigated only two other /3 values, 6.0 and 6.4. A more complete study of SU(3) 
data is in preparation [25], and will provide a more stringent test of scaling. 

When comparing quenched and unqucnched QCD we clearly see that the effect 
of adding fermions can not be absorbed into a "/3-shift". If one matches the 
string tensions in this manner, one finds that the running coupling in the UV is 
quite different (sect. 4.1.1). Also, our results are accurate enough to show, for 
the first time, that the dimensionless quantity ro^/a is different in quenched and 
unquenched QCD (sect. 5.3). 

Till now wc performed only uncorrclatcd fits, ignoring the correlations between 
the potential at different r- values. Except for the case of the SU(2) potential [6], 
the data we used are not strongly correlated. In such cases previous experience 
with (modified) Coulomb + linear fits [23, 27, 28] indicates, that neither the fit 
parameters nor their errors change significantly when employing correlated fits. We 
have explicitly checked this (sect. 4.1.2) by performing fits of the potential data to 
modified Coulomb + linear ansatze and comparing them with previous results from 
correlated fits. For strongly correlated data we find that our naive uncorrelated 
fit over- not underestimates the errors (what uncorrelated fits do underestimate, 
is x^)- If is easy to understand why this is true. So in all cases our errors seem to 
be on the safe side. 

However, what clearly should be done, is to merge our fitting routine with that 
for the potential and perform fits on bootstrap copies of the potential; this would 
automatically take into account all correlations in an optimal manner and simplify 
the error analysis. Note that once this has been done, another nice feature of 
our approach can be fully exploited. Namely, that our method just needs one fit 
of potential data to obtain the running coupling aylq) at all (accessible) scales; 
in particular, we obtain estimates of the UV coupling av{q*) (or Ay, if desired), 
the string tension a, and ro in one fell swoop. At present one usually performs a 
separate fit (or even simulation) for each of these quantities. 

We know how to incorporate QCD to any number of loops in our ansatz for the 
running coupling (sect. 2.1). Once the third /3-function coefficient, P2, becomes 
a\'ailal)le in the V scheme it should of course be incorporated into our ansatz. 

■^^We carefully checked that the one potential problem with global fits, systematic errors from 
the form of the ansatz, does not lead to significant errors in the intermediate region where ro is 
situated. 
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With present data we can not determine this parameter rehably by fitting. It is 

suggestive, though, that pure SU(3) on finer lattices seems to favor a value of Oi 
(which is related to j32) that is in the neighborhood of the value this quantity has 
in the MS scheme (cf. sect. 4.1). With future high precision data it therefore might 
be possible to determine P2 by fitting. On the other hand, non-perturbative effects 
might become important before 3-loop effects are clearly visible. 

Besides fitting lattice potential data and extracting the running coupling along 
the way, our results should have various other applications in the upcoming era 
of high precision Monte Carlo calculations. For example, Monte Carlo simulations 
are now competing (see [35] for a recent review) with — or have already surpassed 
— experiment in determining the running coupling of real world QCD at high 
energy (this is due to the fact that in the simulation one can use a more suitable 
observable, like level sphttings in bottonium [34, 35] , to set the scale than is possible 
in accelerator experiments). To minimize the systematic error in these calculations 
it is crucial to be able to extrapolate the running coupling to different energy scales 
and to a number of fcrmions different from those used in a specific simulation. 

Another obvious application is to potential models (on the lattice or in the 
continuum). As far as we know, previously no simple analytic expression incorpo- 
rating both 2-Ioop QCD and the leading and subleading string pictiuc prediction 
for the potential was known. Although the relation between the static potential 
and heavy quark potentials is not well understood, in high precision studies of the 
heavy quark spectrum our form of the potential might also present an improve- 
ment. More generally, our parameterization of the running coupling should be 
useful for just about any phenomena where one would like to extend perturbation 
theory without running into the Landau pole. 

Finally, we note that our scheme can accommodate string breaking in the pres- 
ence of light quarks. If the onset of string breaking (before the Yukawa-like poten- 
tial between the pair-created mesons becomes relevant) is modeled by a potential 
V{r) = ^(1 — e"'*'') for large r, the small q behavior of the running coupling is 



To take this into account we add a suitable term to eq. (18), 



+ ^^y^ + ln[l + ^ In (co + ^A(, 

(41) 

where there is some freedom in the exact choice of the second term on the rhs. 

Another, less trivial extension would be to generalize our ansatz for the infinite 
lattice potential to one for a finite lattice. It is well known, at least empirically. 
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that the right potential on a finite lattice is not obtained from some infinite lattice 
expression by simply replacing the integral over the Brillouin zone by a corre- 
sponding sum over the lattice. In understanding this problem it will be necessary 
to take into account the exact definition of the finite-volume potential, e.g. via 
Wilson loops, a consideration that played no role in the present work. 
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